Cargamos los datos de ejemplo children de la librería FactoMineR. Seleccionamos solamente las columnas relacionadas con el nivel de estudios de los encuestados.
data("children")
children = children[,1:5]
children
## unqualified cep bepc high_school_diploma university
## money 51 64 32 29 17
## future 53 90 78 75 22
## unemployment 71 111 50 40 11
## circumstances 1 7 5 5 4
## hard 7 11 4 3 2
## economic 7 13 12 11 11
## egoism 21 37 14 26 9
## employment 12 35 19 6 7
## finances 10 7 7 3 1
## war 4 7 7 6 2
## housing 8 22 7 10 5
## fear 25 45 38 38 13
## health 18 27 20 19 9
## work 35 61 29 14 12
## comfort 2 4 3 1 4
## disagreement 2 8 2 5 2
## world 1 5 4 6 3
## to_live 3 3 1 3 4
Creamos a continuación las matrices de frecuencias relativas, condicionadas a filas, a columnas, y marginales de filas y columnas:
miF = children/sum(children); round(miF, 4)
## unqualified cep bepc high_school_diploma university
## money 0.0308 0.0386 0.0193 0.0175 0.0103
## future 0.0320 0.0543 0.0470 0.0452 0.0133
## unemployment 0.0428 0.0669 0.0302 0.0241 0.0066
## circumstances 0.0006 0.0042 0.0030 0.0030 0.0024
## hard 0.0042 0.0066 0.0024 0.0018 0.0012
## economic 0.0042 0.0078 0.0072 0.0066 0.0066
## egoism 0.0127 0.0223 0.0084 0.0157 0.0054
## employment 0.0072 0.0211 0.0115 0.0036 0.0042
## finances 0.0060 0.0042 0.0042 0.0018 0.0006
## war 0.0024 0.0042 0.0042 0.0036 0.0012
## housing 0.0048 0.0133 0.0042 0.0060 0.0030
## fear 0.0151 0.0271 0.0229 0.0229 0.0078
## health 0.0109 0.0163 0.0121 0.0115 0.0054
## work 0.0211 0.0368 0.0175 0.0084 0.0072
## comfort 0.0012 0.0024 0.0018 0.0006 0.0024
## disagreement 0.0012 0.0048 0.0012 0.0030 0.0012
## world 0.0006 0.0030 0.0024 0.0036 0.0018
## to_live 0.0018 0.0018 0.0006 0.0018 0.0024
margFilas = rowSums(miF); margFilas
## money future unemployment circumstances hard
## 0.116405308 0.191797346 0.170687575 0.013268999 0.016284680
## economic egoism employment finances war
## 0.032569361 0.064535585 0.047647768 0.016887817 0.015681544
## housing fear health work comfort
## 0.031363088 0.095898673 0.056091677 0.091073583 0.008443908
## disagreement world to_live
## 0.011459590 0.011459590 0.008443908
margCols = colSums(miF); margCols
## unqualified cep bepc high_school_diploma
## 0.19963812 0.33594692 0.20024125 0.18094089
## university
## 0.08323281
condFilas = miF/margFilas; round(condFilas, 4); rowSums(condFilas) # Matriz R
## unqualified cep bepc high_school_diploma university
## money 0.2642 0.3316 0.1658 0.1503 0.0881
## future 0.1667 0.2830 0.2453 0.2358 0.0692
## unemployment 0.2509 0.3922 0.1767 0.1413 0.0389
## circumstances 0.0455 0.3182 0.2273 0.2273 0.1818
## hard 0.2593 0.4074 0.1481 0.1111 0.0741
## economic 0.1296 0.2407 0.2222 0.2037 0.2037
## egoism 0.1963 0.3458 0.1308 0.2430 0.0841
## employment 0.1519 0.4430 0.2405 0.0759 0.0886
## finances 0.3571 0.2500 0.2500 0.1071 0.0357
## war 0.1538 0.2692 0.2692 0.2308 0.0769
## housing 0.1538 0.4231 0.1346 0.1923 0.0962
## fear 0.1572 0.2830 0.2390 0.2390 0.0818
## health 0.1935 0.2903 0.2151 0.2043 0.0968
## work 0.2318 0.4040 0.1921 0.0927 0.0795
## comfort 0.1429 0.2857 0.2143 0.0714 0.2857
## disagreement 0.1053 0.4211 0.1053 0.2632 0.1053
## world 0.0526 0.2632 0.2105 0.3158 0.1579
## to_live 0.2143 0.2143 0.0714 0.2143 0.2857
## money future unemployment circumstances hard
## 1 1 1 1 1
## economic egoism employment finances war
## 1 1 1 1 1
## housing fear health work comfort
## 1 1 1 1 1
## disagreement world to_live
## 1 1 1
condCols = t(t(miF)/margCols); round(condCols, 4); colSums(condCols) # Matriz C
## unqualified cep bepc high_school_diploma university
## money 0.1541 0.1149 0.0964 0.0967 0.1232
## future 0.1601 0.1616 0.2349 0.2500 0.1594
## unemployment 0.2145 0.1993 0.1506 0.1333 0.0797
## circumstances 0.0030 0.0126 0.0151 0.0167 0.0290
## hard 0.0211 0.0197 0.0120 0.0100 0.0145
## economic 0.0211 0.0233 0.0361 0.0367 0.0797
## egoism 0.0634 0.0664 0.0422 0.0867 0.0652
## employment 0.0363 0.0628 0.0572 0.0200 0.0507
## finances 0.0302 0.0126 0.0211 0.0100 0.0072
## war 0.0121 0.0126 0.0211 0.0200 0.0145
## housing 0.0242 0.0395 0.0211 0.0333 0.0362
## fear 0.0755 0.0808 0.1145 0.1267 0.0942
## health 0.0544 0.0485 0.0602 0.0633 0.0652
## work 0.1057 0.1095 0.0873 0.0467 0.0870
## comfort 0.0060 0.0072 0.0090 0.0033 0.0290
## disagreement 0.0060 0.0144 0.0060 0.0167 0.0145
## world 0.0030 0.0090 0.0120 0.0200 0.0217
## to_live 0.0091 0.0054 0.0030 0.0100 0.0290
## unqualified cep bepc high_school_diploma
## 1 1 1 1
## university
## 1
Nota: A partir de las matrices R o C, podríamos ver si hay filas o columnas con perfiles similares y agruparlas según el principio de equivalencia distribucional.
EJERCICIO 1
Crea un gráfico que permita visualizar las similitudes entre los perfiles dados por las frecuencias condicionadas a las filas (matriz R) y crea una nueva tabla de contingencia agrupando las dos categorías más similares.
Aplicamos el test de independencia de la para contrastar la hipótesis nula de independecia entre las variables “nivel de estudios” y “razones para no tener hijos”.
chisq.test(children)
## Warning in chisq.test(children): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test
##
## data: children
## X-squared = 123.82, df = 68, p-value = 4.145e-05
El resultado del test indica que existe una clara dependencia entre las variables estudiadas. Se genera un warning porque algunas de las frecuencias absolutas son inferiores a 5, que es lo recomendado para aplicar este test. No obstante, a la vista del p-valor, podemos asumir la dependencia. Estudiaremos la naturaleza y causas de dicha dependencia con un AFC simple.
Generaremos un modelo AFC con la librería FactoMineR y seleccionaremos el número óptimo de dimensiones (o componentes) del modelo. Con el número de dimensiones elegido, estimaremos el modelo AFC final.
res.afc = CA(children, graph = FALSE)
eig.val <- get_eigenvalue(res.afc)
Vmedia = 100 * (1/nrow(eig.val))
fviz_eig(res.afc, addlabels = TRUE) +
geom_hline(yintercept=Vmedia, linetype=2, color="red")
kable(eig.val)
| eigenvalue | variance.percent | cumulative.variance.percent | |
|---|---|---|---|
| Dim.1 | 0.0387565 | 51.89444 | 51.89444 |
| Dim.2 | 0.0188857 | 25.28778 | 77.18222 |
| Dim.3 | 0.0088874 | 11.90017 | 89.08240 |
| Dim.4 | 0.0081536 | 10.91760 | 100.00000 |
res.afc = CA(children, graph = FALSE, ncp = 2)
Seleccionaremos 2 componentes (dimensiones), que explican un 77.2% del total de variabilidad de los datos.
EJERCICIO 2
¿Cuál es el número máximo de dimensiones que se podrían obtener en este estudio? ¿Por qué?
A continuación, generaremos algunos ejemplos de gráficos que se pueden utilizar para interpretar los resultados a nivel de filas del modelo AFC, es decir, para ver las categorías (en filas) que son similares o diferentes entre sí, y cuáles han contribuido más a crear cada dimensión.
head(res.afc$row$contrib) # Contribución relativa de cada fila a cada dimensión
## Dim 1 Dim 2
## money 4.142727 2.715751
## future 9.791562 21.646842
## unemployment 23.965450 2.029596
## circumstances 6.237428 2.944103
## hard 2.529531 1.002058
## economic 12.539185 12.431988
head(res.afc$row$cos2) # Calidad de la representación de cada fila en cada dimensión
## Dim 1 Dim 2
## money 0.4260428 0.13609637
## future 0.4606709 0.49627585
## unemployment 0.9436294 0.03894173
## circumstances 0.7220943 0.16608530
## hard 0.8079044 0.15595614
## economic 0.6454456 0.31183200
fviz_ca_row(res.afc, axes = c(1,2), repel = TRUE, col.row = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
fviz_contrib(res.afc, choice = "row", axes = 1)
fviz_contrib(res.afc, choice = "row", axes = 2)
fviz_cos2(res.afc, choice = "row", axes = 1)
fviz_cos2(res.afc, choice = "row", axes = 2)
fviz_cos2(res.afc, choice = "row", axes = 1:2)
corrplot(res.afc$row$cos2[,1:2], is.corr=FALSE)
EJERCICIO 3
¿Qué significado le podríamos dar a cada dimensión según la variable en las filas? ¿Qué información de dicha variable resume cada dimensión?
Generaremos también algunos de los gráficos que se pueden utilizar para interpretar los resultados a nivel de columnas del modelo AFC.
head(res.afc$col$contrib) # Contribución relativa de cada columna a cada dimensión
## Dim 1 Dim 2
## unqualified 26.439479 0.02620065
## cep 15.522144 3.40369039
## bepc 3.203853 7.40770166
## high_school_diploma 31.380527 24.39373950
## university 23.453997 64.76866780
head(res.afc$col$cos2) # Calidad de la representación de cada columna en cada dimensión
## Dim 1 Dim 2
## unqualified 0.6934614 0.0003348661
## cep 0.5701157 0.0609187629
## bepc 0.1429008 0.1610034484
## high_school_diploma 0.6382066 0.2417515320
## university 0.4206562 0.5660635034
fviz_ca_col(res.afc, axes = c(1,2), repel = TRUE, col.col = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
fviz_contrib(res.afc, choice = "col", axes = 1)
fviz_contrib(res.afc, choice = "col", axes = 2)
fviz_cos2(res.afc, choice = "col", axes = 1:2)
corrplot(res.afc$col$cos2[,1:2], is.corr=FALSE)
EJERCICIO 4
¿Qué significado le podríamos dar a cada dimensión según la variable en las columnas? ¿Qué información de dicha variable resume cada dimensión?
El principal objetivo del AFC era entender la relación entre las dos variables estudiadas. Para ello, podemos comparar los gráficos de filas y columnas obtenidos anteriormente o bien utilizar un biplot, es decir, un gráfico que nos permite representar conjuntamente las filas y columnas de la tabla de contingencia.
Existen diferentes formas de representar un biplot en el AFC simple. En un biplot simétrico, solo tiene sentido interpretar la distancia entre filas o la distancia entre columnas pero no la distancia entre una fila y una columna. Solo se pueden hacer afirmaciones generales del patrón de la relación entre filas y columnas. Por tanto, si queremos estudiar las relaciones entre filas y columnas, será más adecuado un biplot asimétrico. Existen muchos tipos de biplots asimétricos, representaremos a continuación uno de ellos.
# Simetricos
fviz_ca_biplot(res.afc, repel = TRUE)
# Asimetricos
fviz_ca_biplot(res.afc,
map ="rowprincipal", repel = TRUE) # Columnas representadas en espacio de las filas
## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
fviz_ca_biplot(res.afc,
map ="colprincipal", repel = TRUE) # Filas representadas en espacio de las columnas
EJERCICIO 5
A la vista de todos los resultados y gráficos generados, ¿qué puedes afirmar sobre la dependencia de las variables estudiadas? ¿a qué se debe principalmente que exista esa dependencia? ¿opinan de igual modo sobre las razones de no tener hijos los encuestados con cualquier nivel de estudios?
Cargamos la base de datos. Descartaremos la columna Time, puesto que desconocemos su significado exacto y su relevancia en el análisis. En primer lugar, generaremos un gráfico para visualizar el número de observaciones en cada categoría de las variables a incluir en el modelo AFC, por si hubiera alguna variable constante o casi constante que debamos eliminar.
data(poison)
head(poison)
## Age Time Sick Sex Nausea Vomiting Abdominals Fever Diarrhae Potato
## 1 9 22 Sick_y F Nausea_y Vomit_n Abdo_y Fever_y Diarrhea_y Potato_y
## 2 5 0 Sick_n F Nausea_n Vomit_n Abdo_n Fever_n Diarrhea_n Potato_y
## 3 6 16 Sick_y F Nausea_n Vomit_y Abdo_y Fever_y Diarrhea_y Potato_y
## 4 9 0 Sick_n F Nausea_n Vomit_n Abdo_n Fever_n Diarrhea_n Potato_y
## 5 7 14 Sick_y M Nausea_n Vomit_y Abdo_y Fever_y Diarrhea_y Potato_y
## 6 72 9 Sick_y M Nausea_n Vomit_n Abdo_y Fever_y Diarrhea_y Potato_y
## Fish Mayo Courgette Cheese Icecream
## 1 Fish_y Mayo_y Courg_y Cheese_y Icecream_y
## 2 Fish_y Mayo_y Courg_y Cheese_n Icecream_y
## 3 Fish_y Mayo_y Courg_y Cheese_y Icecream_y
## 4 Fish_y Mayo_n Courg_y Cheese_y Icecream_y
## 5 Fish_y Mayo_y Courg_y Cheese_y Icecream_y
## 6 Fish_n Mayo_y Courg_y Cheese_y Icecream_y
poison2 <- poison[, -2]
# Gráfico con las frecuencias de cada categoría para todas las variables
frecus = apply(poison2[,4:14], 2, table)
rownames(frecus) = c("No", "Yes")
barplot(frecus, las = 2, legend = TRUE, cex.names = 0.8)
Deberíamos plantearnos eliminar la variable Fish, porque casi todos los individuos han tomado pescado. La mantendremos en la base de datos porque el pescado puede ser una probable causa de intoxicación alimentaria y no nos interesa perder esta información.
A continuación, vamos a realizar un análisis descriptivo sencillo de las variables auxiliares:
summary(poison2$Age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 6.00 8.00 16.93 10.00 88.00
barplot(table(poison2$Age), cex.names = 0.7)
table(poison2$Sick)
##
## Sick_n Sick_y
## 17 38
table(poison2$Sex)
##
## F M
## 28 27
Detectamos anomalías en la variable edad, ya que hay edades desde 36 a 88 años y se trataba de niños de primaria. Vamos a suponer que estas edades se han introducido mal y que realmente les falta el símbolo decimal entre los dos dígitos, por lo que corregimos estos valores.
poison2$Age[poison2$Age >= 36] = poison2$Age[poison2$Age >= 36]/10
hist(poison2$Age, col = "grey", xlab = "Edad", main = "Histograma")
Generaremos un modelo AFC múltiple, seleccionaremos el número de componentes óptimo y generaremos de nuevo el modelo para ese número de dimensiones óptimo.
res.mca <- MCA(poison2, graph = FALSE, quanti.sup = 1, quali.sup = 2:3)
eig.val <- get_eigenvalue(res.mca)
Vmedia = 100 * (1/nrow(eig.val))
fviz_eig(res.mca, addlabels = TRUE) +
geom_hline(yintercept=Vmedia, linetype=2, color="red")
kable(head(eig.val))
| eigenvalue | variance.percent | cumulative.variance.percent | |
|---|---|---|---|
| Dim.1 | 0.3352314 | 33.523140 | 33.52314 |
| Dim.2 | 0.1291398 | 12.913979 | 46.43712 |
| Dim.3 | 0.1073485 | 10.734849 | 57.17197 |
| Dim.4 | 0.0958795 | 9.587950 | 66.75992 |
| Dim.5 | 0.0788328 | 7.883277 | 74.64319 |
| Dim.6 | 0.0710898 | 7.108981 | 81.75217 |
res.mca <- MCA(poison2, graph = FALSE, quanti.sup = 1, quali.sup = 2:3, ncp = 3)
Seleccionaremos 3 dimensiones, que explican el 57.2% de la inercia total.
A continuación se muestran algunos gráficos ejemplo para poder interpretar el modelo. Es importante recordar que es necesario explorar TODAS las dimensiones seleccionadas, aunque aquí por brevedad no se repitan siempre los gráficos para las 3 dimensiones. Tampoco significa que SIEMPRE tengamos que hacer todos los tipos gráficos. En cada estudio deberemos elegir los gráficos más representativos e informativos, la forma de colorearlos o mostrarlos, etc.
Es importante también que se mantenga la misma escala en los dos ejes representados.
# Podemos colorear los individuos segun los valores de cualquier
# variable (activa o suplementaria) en nuestros datos, por ejemplo Sick
fviz_mca_ind(res.mca, axes = c(1,2),
label = "none", # ocultar etiquetas de los individuos
habillage = "Sick", # variable utilizada para colorear
palette = c("#00AFBB", "#E7B800"), # colores para cada grupo
addEllipses = TRUE, # elipse alrededor de cada grupo coloreado
ellipse.type = "confidence", # elipse de confianza alrededor del punto medio
ggtheme = theme_minimal())
# Coloreamos por una variable numérica como Age, ahora representando las dimensiones 1 y 3
fviz_mca_ind(res.mca, axes = c(1,3),
label = "none",
col.ind = poison2$Age, # variable utilizada para colorear
ggtheme = theme_minimal())
# Para colorear individuos según varias variables categóricas
fviz_ellipses(res.mca, c("Fever", "Sick"), geom = "point")
# Correlación entre variables y dos primeras dimensiones
fviz_mca_var(res.mca, choice = "mca.cor",
repel = TRUE,
ggtheme = theme_minimal())
# Correlación entre variables y dimensiones
corrplot(res.mca$var$cos2, is.corr=FALSE, tl.col = 1, win.asp = 0.5, cl.ratio = 0.3, cl.cex = 0.8)
fviz_cos2(res.mca, choice = "var", axes = 1:3)
fviz_mca_var(res.mca, col.var="steelblue", shape.var = 15, repel = TRUE)
## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
# Contributions de las variables a la dimension 1
fviz_contrib(res.mca, choice = "var", axes = 1, top = 15)
# Contributions de las variables a la dimension 2
fviz_contrib(res.mca, choice = "var", axes = 2, top = 15)
# Contributions de las variables a la dimension 3
fviz_contrib(res.mca, choice = "var", axes = 3, top = 15)
# Variables con cos2 >= 0.4
fviz_mca_var(res.mca, select.var = list(cos2 = 0.4))
# Top 10 variables con el cos2 mas alto
fviz_mca_var(res.mca, select.var= list(cos2 = 10))
# Incluir en el gráfico solo un conjunto de variables seleccionadas
misvar <- list(name = c("Fever_n", "Abdo_y", "Diarrhea_n",
"Fever_Y", "Vomit_y", "Vomit_n"))
fviz_mca_var(res.mca, select.var = misvar)
fviz_mca_biplot(res.mca, geom.ind = "point",
repel = TRUE, # Avoid text overlapping (slow if many point)
ggtheme = theme_minimal())
## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
# top 10 individuos y variables con contribuciones mas altas
fviz_mca_biplot(res.mca, select.ind = list(contrib = 10),
select.var = list(contrib = 10), repel = TRUE,
ggtheme = theme_minimal())
EJERCICIO 6
A la vista de todos los resultados y gráficos anteriores, ¿qué tipo de información resume la primera y la segunda componente? ¿qué relación existe entre los alimentos ingeridos y los síntomas presentados?